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ABSTRACT 


A rapid implicit-explicit finite difference scheme has been developed 
for the numerical solution of the time-dependent, incompressible, two- 
dimensional Navier-Stokes equations in conservation-law form using vorticity 
and stream function variables. The algorithm is second-order time-accurate 
and spatially factored. The systems of equations are solved at each time 
step by an iterative technique. Numerical results have been obtained with 
the technique for a circular cylinder at a Reynolds number of 15, and an 
NACA 0012 airfoil at zero angle of attack at Reynolds numbers of 10 3 and 10 4 . 
The results are in agreement with another numerical technique, and the com- 
puting time required to obtain the steady-state solution at the Reynolds 
number of 10 4 was 49.7 sec on a CDC 7600 computer using a 65 x 84 computa- 
tional grid. 


INTRODUCTION 


Rapid numerical solutions to the time-dependent Navier-Stokes equations 
at moderate Reynolds numbers and low Mach numbers have been difficult to 
obtain because of the elliptic nature of the equations in this regime. 

Unlike high Mach number compressible flows, each portion of the flow field 
significantly influences every other portion of the flow field which dras- 
tically increases the amount of computation time required to obtain either 
a time-accurate or steady-state numerical solution. 

Fast numerical techniques developed for the solution of the compressible 
Navier-Stokes equations, such as the algorithms of MacCormack 1 and Beam and 
Warming, 2 unfortunately do not perform as efficiently in the low Mach number 
regime as they do at high Mach numbers. J. L. Steger (private communication) 
has applied the latter technique to the laminar flow (Re = 10 4 ) about an 
NACA 0012 airfoil at zero angle of attack and a Mach number of 0.2 and was 
able to obtain a steady-state solution on a 71 x 33 grid in about 20 min of 
computing time on a CDC 7600 computer. 

Roache 3 has presented a fast numerical technique for the incompressible 
steady Navier-Stokes equations as applied to the internal flow in flush 
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inlets. However, his program at present does not simulate unsteady flows 
and has only been tested on a very coarse mesh, where the number of itera- 
tions required to converge to a steady-state solution may be lower than the 
number required for a fine mesh. Other excellent numerical techniques for 
solving the incompressible time-dependent Navier-Stokes equations, such as 
the methods of Mehta, 1 ® Reddy and Thompson, 1 * Wu, 5 and Thames, 6 all require 
approximately one-half hour or more of computing time on a CDC 7600 to obtain 
steady-state solutions for moderate Reynolds numbers (-10 1 *) and sufficiently 
fine grids (-5000 mesh points). Even the very efficient method of Ghia, 
Hankey, and Hodge, 7 which was developed for cavity flows, appears to require 
about 4 min of computing time at these Reynolds numbers and grid sizes to 
obtain steady-state solutions. It is therefore the purpose of this report 
to present a rapid numerical technique for solving the incompressible time- 
dependent Navier-Stokes equations for the flow about two-dimensional airfoils 
using very fine computational grids. 

The present numerical technique evolved as the result of an investiga- 
tion intended to obtain numerical solutions to the incompressible time- 
dependent Navier-Stokes equations by utilizing the Newton-Kantorovich tech- 
nique . 9 This technique is a rigorous mathematical concept based on functional 
analysis that transforms a nonlinear partial differential equation into a 
sequence of linear partial differential equations. When this technique is 
applied to a second-order time-accurate differencing of the Navier-Stokes 
equations and maximum advantage is taken of previous time step information, 
it reduces to approximating the Navier-Stokes equations by a single linear 
partial differential equation. The equation becomes linear because the non- 
linear term has, in effect, been replaced by a first-order Taylor series 
expansion in function space. A detailed discussion of the mathematical 
foundations of the technique as well as a summary of previous investigations 
utilizing the technique are given in Gabrielsen 8 and Davis, Gabrielsen, and 
Mehta . 9 

As a result of that investigation, new insight was gained to improve the 
numerical technique used to solve the Navier-Stokes equations. This report 
describes these improvements and the results obtained to date using the 
improved technique. 


MATHEMATICAL FORMULATION 


The present numerical technique has been incorporated into a computer 
code that is restricted to symmetric flow fields and geometries. It is a 
special version of a more general computer code that was developed for non- 
symmetric airfoils at angle of attack and is described in detail in refer- 
ence 10. The mathematical formulation contained herein describes only the 
symmetric case, although the development for the angle of attack case is 
very similar. 

The flow field and geometry under consideration are mapped into the unit 
half-circle. As a result of this transformation, the equation governing the 
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unsteady, incompressible flow of a Newtonian fluid may be expressed in terms 
of a transformed vorticity transport equation and disturbance stream function 
equation as 
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and 

w = vorticity 

y = total stream function 

^ = disturbance stream function 'H-y 

(x,y) = coordinates of physical flow field 

(r,0) = coordinates of the interior of the unit half-circle where 0 < r < 1 
and tt < 0 < 2tt 

p = stretched radial coordinate where the stretching function is defined 

below 

L = dimensionless chord 

R = Reynolds number Ul/v 

U = freestream velocity 

1 = chord 

v = kinematic viscosity 

t = time 
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is the conservation-law form of the convective terms. 
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The radial coordinate r is transformed according to a function that 
stretches the region near the body surface (r = 1) so that proper resolution 
of the boundary layer and separation zone may be obtained. The stretched 
radial coordinate is depicted by the variable p and is expressed as 


P = (k 2 + k 2 ) 1 [tanh -1 (rk 3 - k 4 ) + k 2 ] (5) 

with 

k 3 = jtanh(k x ) + tanh(k 2 )J £l - — ^T)"] 


and 


k 4 = tanh(k 2 ) - r Q [tanh(k x ) + tanh(k 2 )] 



1 ) 


The constants r Q , k x , and k 2 (all positive) determine the value of p. As 
r varies from r G to 1, p varies from 0 to 1. 


The mapping into the unit circle is accomplished by the transformation 
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A proper choice of the constants y and c invokes the solution for the flow 
over any one of a number of shapes, including a flat plate, a circular or 
elliptical cylinder, or thick airfoils. The trailing edge of an airfoil 
shape may be rounded off by defining 


c = |? + (1 - n 2 ) 1/,2 J(l - <$) , where 0 < 6 < 1 (7) 

The components of velocity u x and u 2 are defined in terms of the dis- 
turbance stream function as 


U 1 


- 1 _ (M. , 1 I_\ 
rH \30 30/ ’ 


u 2 


1. /dp_ djh 
H Vdr 3p + 



( 8 ) 


For Reynolds numbers much larger than unity, the vorticity in the flow 
field exists only near the body and in the wake. Away from this region, the 
flow is essentially irrotational; therefore, the region of calculation is 
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subdivided into two parts: a small viscous region and a large irrotational 

region bounded by and p^_ with tt < 0 < 0 2 (see fig. 1) . 

The boundary conditions applied along the boundaries of the region of 
calculation are shown in figure 1 as follows: 

1. On the surface (p = 1) , the constraint of no slip is applied in the 

form 


T = 0 (or ip = -y) 


(9) 


and 



( 10 ) 


Condition (10) is used to calculate the surface vorticity by the method of 
Israeli-'- described in detail later. 


2. Along the line of symmetry, the vorticity w and the disturbance 
stream function if' are specified to be zero. 

3. The flow at the far boundary is constrained with first-order differ- 
ential relations obtained from the Navier-Stokes equations by dropping the 
tangential derivative of the pressure and viscous terms (i.e., at the outer 
boundary); the inertia terms are dominant: 
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and u 2 is obtained from the 0 component of the Navier-Stokes equation. 
That is , 
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Note that, in equation (13), w = 0 when it < 0 < 0 2 . At t = 0, the flow is 
irrotational (without circulation), that is, m = 0 and ip = -yr 2 . 


These boundary conditions are believed to be superior to specifying 
either potential flow or uniform velocity since eddies or vortices can pass 
through the downstream boundary. Also, because the velocity far away is not 
defined, the circulation there can change with time. Therefore, equation (11) 
correctly represents the vorticity transport through the downstream boundary. 
In equation (13), the absence of the tangential pressure derivative will not 
significantly affect the motion of a vortex through the boundary. 
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The surface pressure distribution is obtained by integrating the tangen- 
tial component of the Navier-Stokes equation. That is, 

e-Mf r$ dS (14) 

O K 

where p(0) = 0. The pressure coefficient Cp is, therefore, equal to 2p. 

On the surface, the tangential stress is given by a\2 = (L/R)m. Both p and 
a 12 are made dimensionless with pU 2 . The coefficients of lift and moment 
are zero for symmetric geometries and flow fields; however, the coefficients 
of drag are computed as 
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where the subscripts P and F, respectively, represent the contributions due 
to pressure and viscous forces. 


Equations (1) through (4) represent a nonlinear coupled set of partial 
differential equations. To reduce the difficulty in solving these equations, 
the nonlinearity may be removed by applying the Newton-Kantorovich transfor- 
mation to the equations 9 or by performing a Taylor series expansion in func- 
tion space to the nonlinear terms. Employing such a procedure alters the 
convection terms in the vorticity transport equation, and the linear approxi- 
mation to this equation becomes 
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where 


and 


A to = a) — w 
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(19) 
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The functions w 0 and Y 0 are initial guesses to the solutions us and Y . 
These initial guesses must be sufficiently accurate such that the truncation 
error associated with approximating the nonlinear terms is no worse than the 
truncation error of the overall numerical technique. 


NUMERICAL FORMULATION 


The Navier-Stokes equations were solved by a three-point backward time 
differencing and implicit-explicit factored central space differencing scheme. 
The spatially factored, time— differenced expression for the linearized approxi- 
mation to the vorticity transport equation^ is 
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and the cross-derivative terms Q(aj n ,Y ), F(co ,Y ,Y n ) , and G(to ,Y ,Y n ) are 
given as 
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The values of the three-point backward time difference parameters are T = 3, 
= -4, and T 2 = 1. The transformation parameter A is 


A = H 2 r 2 I 1 
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Let the 

initial guesses at the solution m and T be given by 
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at the first two time steps where they are given as 
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except at the first two time steps. Therefore, except at the first two time 
steps, the terms F and G may be neglected because they are the same order 
as the overall truncation error. In addition, the second-order term on the 
right-hand side of equation (20) is made explicit by approximating it by 
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The final second-order accurate time-differenced expression for the vorticity 
transport equation is then 
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Consistent with the spatial factoring concept, equation (34) is split into 
two equations whose finite difference analogies each produce tridiagonal sys- 
tems of equations. That is. 
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The spatial derivatives associated with equations (35) and (36) , as well as 
equation (19) , were approximated by central differencing formulas everywhere 
except at the boundaries of the flow field. The truncation error for the 
vorticity transport equation is 0[(Ap) 2 + (A0) 2 + (At) 2 ], except at the 
first two time steps where the temporal error is O(At) , because a two-point 
backward difference formula is required such that T = 2, Tj = -2, and 
T 2 = 0. 


The disturbance stream function equation is solved by the method of ref- 
erence 12. The details of the procedure are not described here, but the 
truncation error of the disturbance stream function equation is 
0[ ( Ap) 2 + (A6) 2 ] . 


At the outer boundary, the finite difference forms of equations (35) 
and (36) are used to describe the vorticity condition, except that the diffu- 
sion terms are deleted and forward space differencing is used in the p 
direction. This formulation results in a first-order truncation error in the 
space variable p. The disturbance stream function at this boundary is 
obtained from the finite-difference form of equations (12) and (13) : 



1 

3 





+ 4ip 


i,2 



(37) 




9 



where 
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Along the lines of symmetry the boundary conditions w = 0 and ip = 0 are 
specified. The no-slip boundary condition represented by equations (9) 
and (10) is satisfied by specifying the condition ip = 0 along the boundary 
and by iteratively determining the surface vorticity that numerically 
enforces dip / 3 p = 0. The latter condition is accomplished by employing the 
method of Israeli. 11 That is, the surface vorticity at a particular itera- 
tion "k" is related to dip/ 9p and the surface vorticity at the previous 
iteration "k-1" by 
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The coefficient g is a relaxation parameter whose optimum value is deter- 
mined in reference 11 for boundary-layer type flows. For the fine grids and 
relatively large time steps used in the present study, the equations defining 
the optimum value of g as given in reference 11 may be approximated to the 

Ar jL 

first order in as 

yvAt 


g 
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where Atjl is the normal distance from the surface to the next grid point 
and v, R, and At are the kinematic viscosity, Reynolds number, and time 
step size, respectively. The value of g for a given geometry and computa- 
tional mesh must be determined manually once for a particular Reynolds num- 
ber and time step size. Then, to determine the value of g at any other 
Reynolds number or time step size, equation (40) may be used with reasonable 
accuracy. In fact, figure 2 is a log-log plot of g versus At for a NACA 
0012 airfoil at a Reynolds number of 10 4 using a 65 x 84 computational mesh. 
Equation (40) is represented in the figure by the straight line and optimum 
values of g determined manually are shown by the symbols. Note that the 
agreement is reasonable. In fact, the same calculations using the manually 
determined g required only a few percent less compute time than the value 
of g determined by equation (40) . Furthermore, the savings in engineering 
and compute time by not manually searching for the optimum value of g more 
than makes up for the difference. More accurate formulations of the optimum 
value of g could be determined. However, for the purposes of the current 
study, this did not appear to be warranted. 
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In solving the systems of equations described by equations (2), (35), 
and (36), together with the boundary conditions (37), (38), and (39), a mixed 
iterative/noniterative procedure is used. Iteration is required because w n 
is contained in the term Q(( 0 n ,'l' o ) on the right-hand side of equation (35) . 
Hence, equations (2) , (35) , and (36) must be solved at each time step by 
iterative improvement. In addition, the boundary conditions (37) , (38) , 
and (39) are iterative in nature, hence they must also be satisfied by itera- 
tive improvement at each time step. However, the coefficient matrices of 
equations (35) and (36) do not need to be iteratively improved at each time 
step. A considerable amount of compute time, therefore, is saved at each 
time step by calculating these matrices and performing the forward elimina- 
tion upon them outside the iteration loop. Then, in iteratively improving 
the values of oi n in Q(o3 n ,’i' 0 ) and solving for w* and a new u) n , only back 
substitutions need be performed. Convergence of the above-described itera- 
tion procedure was considered achieved when the change in surface vorticity 
between successive iterations was less than a prescribed value, typically 
10" 3 or 10" 4 . 

Computations of pressure coefficients on the surface and a determination 
of the loads are made at each time step with a finite-difference integration 
formula derived by combining two four-point expressions (in order to have a 
lower effective truncation error) . The details of the formula may be found 
in reference 13. 

DISCUSSION OF RESULTS 


The present numerical technique was used to perform flow-field calcula- 
tions for four test cases. These test cases included a steady-state solution 
for a circular cylinder at Reynolds number 15, a steady-state solution for 
an NACA 0012 airfoil at zero angle of attack at a Reynolds number of 10 3 , and 
a steady-state solution and a time-accurate solution for the same airfoil at 
a Reynolds number of 10 4 . Results were also obtained for the same airfoil at 
a Reynolds number of 10 5 . The results for that test case, however, indicated 
that -the flow field in the separation zone was unsteady at this Reynolds 
number. At this time it is not known whether this unsteadiness is due to 
numerical problems or whether it is a true flow-field phenomena; therefore, 
the results at this Reynolds number are not presented here. 

A summary of the results of all test cases is presented in table 1. 
Included in the table are the compute times required for a steady-state solu- 
tion on a CDC 7600 computer as well as the number and size of the time steps 
used and the size of the computational grid. All test cases were initiated 
from an impulsive start and computed to a steady-state solution. A detailed 
discussion of the results of each test case is given below. 


Circular Cylinder at Re = 15 

The circular cylinder test case was performed on a 33 x 40 grid in the 
half-plane. Plots of the steady-state surface vorticity and pressure 
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coefficient distributions are shown in figure 3. Comparisons of the results 
are made with Mehta's 9 numerical technique. As seen in the figure the 
results are in excellent agreement. Since only a steady-state solution was 
desired for this test case, no attempt was made to obtain time accuracy. As 
seen in table 1 the total compute time required to obtain the solution was 
2.5 sec. A total of eight time steps were taken and the time step sizes (non- 
dimensionalized as UAt/1) were 2.5 diam each. The cylinder traveled a dis- 
tance of 20 diam during this time. The number of iterations required to 
obtain a converged solution at each time step varied from 21 at the first 
time step to 7 at the last time step. 


NACA 0012 Airfoil at Re = 10 3 

A computational grid in the half-plane of 65 x 48 was used to determine 
the flow field about the NACA 0012 airfoil at a Reynolds number of 10 3 . 

Plots of the steady-state surface vorticity and pressure coefficient distri- 
butions are shown in figure 4. As shown in table 1, a total compute time of 
24.8 sec was required to obtain this solution. The size of the time steps 
was initially 0.269 chord. However, the time step size was doubled after 
10 time steps. A total of 23 time steps was used to obtain the steady-state 
solution, and the airfoil traveled a distance of 9.69 chords during this 
time. The number of iterations required for convergence at each time step 
varied from a maximum of 32 at the first time step to a minimum of 8 at the 
last time step. 


NACA 0012 Airfoil at Re = 10 4 

Two test cases were performed using a 65 x 84 computational mesh for the 
NACA 0012 airfoil at a Reynolds number of 10 4 . The results of both test 
cases are shown in figure 5 as plots of the surface vorticity and pressure 
coefficient distributions. The results of both test cases are also compared 
with computations using Mehta's 9 numerical technique. 

The first test case was performed to determine the steady-state, solution 
only. This test case required 49.7 sec of compute time. Since time accuracy 
was not a consideration, the initial time step size was set at 0.138 chord, 
then doubled every 10 time steps. The steady-state solution was obtained 
after the airfoil had traveled about 9.67 chords, hence a total of 30 time 
steps was needed in obtaining the results. The number of iterations required 
for convergence at each time step varied from 32 at the first couple of time 
steps to 10 at the final time step. 

The second test case at this Reynolds number was performed in order to 
obtain a time-accurate solution until a steady state was achieved. Therefore, 
the calculation was initiated with a smaller time step, 0.00276 chord. How- 
ever, as before, the time step size was doubled every ten time steps until 
the airfoil had traveled about 9.88 chord lengths. The total number of time 
steps required to obtain this solution was 84, while the total compute time 
required was only 88 sec. The number of iterations required for convergence 
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of the surface vorticity varied from a maximum of 23 to a minimum of 6. The 
maximum occurred at the first time step and the minimum occurred at the tenth 
time step. Each time the time step was doubled, the number of iterations 
increased then gradually decreased from time step to time step until the time 
step doubled again. In general, fewer iterations were required for conver- 
gence for smaller time steps than for larger ones. 

As seen in figure 5, both test cases were in excellent agreement with 
each other as well as with the solutions obtained using Mehta 's a numerical 
technique . 


CONCLUSIONS 


Based on the results obtained to date, the present method appears capable 
of obtaining rapid steady-state and time-accurate solutions to the incompres- 
sible time-dependent Navier-Stokes equations. It has been successfully applied 
to a circular cylinder at a Reynolds number of 15 and an NACA 0012 airfoil at 
Reynolds numbers of up to 10 4 using fine computational grids. The numerical 
technique itself may be extended to representations of the Navier-Stokes 
equations in the primitive variables u, v, and p and therefore to three- 
dimensional flows. It may also be extended to the time-averaged Navier-Stokes 
equations with turbulence modeling. 
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FIGURE CAPTIONS 



Figure 1.- Domain of calculation, boundary conditions, and grid notation. 

Figure 2.- Comparisons of the optimum value of g as determined manually and 

by equation 40. 

Figure 3.- Surface vorticity and pressure coefficient distributions on a cir- 
cular cylinder; Re = 15. 

Figure 4.- Surface vorticity and pressure coefficient distributions on a 12% 
thick symmetrical airfoil; a = 0; Re = 10 3 . 

Figure 5.- Surface vorticity and pressure coefficient distributions on a 12% 
thick symmetrical airfoil; a = 0; Re = lO 1 *. 
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TABLE 1. SUMMARY OF COMPUTATION TIMES FOR VARIOUS TEST CASES 

USING PRESENT TECHNIQUE 


TEST CASE 

TOTAL COMPUTE 
TIME (CDC 7600), 
seconds 

TOTAL TIME 
STEPS TO 
STEADY STATE 

TIME STEP 

SIZES 

. . . . . 

GRID 

(0 X r) 

CIRCULAR CYLINDER 
AT Re=15 

2.5 

8 

2.5 

(33 X 40) 

NACA 0012 AIRFOIL 
AT Re=10 3 

24.8 

23 


(65 X 48) 

NACA 0012 AIRFOIL 
AT Re=10 4 

49.7 

30 

0.138 TO 0.553 

(65 X 84) 

88.0 

84 

0.00276 TO 0.707 

(65 X 84) 
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